library(amerifluxr) # Package for querying, downloading, and handling AmeriFlux data and metadata.
library(REddyProc) # Package for processing Eddy Covariance flux data (e.g., gap-filling, partitioning)
library(lubridate) # Provides tools to work with date and time (e.g., parsing, formatting).
##
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
library(tidyverse) # A collection of packages for data manipulation, visualization, and more.
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(bigleaf) # Provides functions to calculate leaf area index and other plant trait estimations.
library(data.table) # Offers fast and memory-efficient manipulation of large datasets.
##
## Attaching package: 'data.table'
##
## The following objects are masked from 'package:dplyr':
##
## between, first, last
##
## The following object is masked from 'package:purrr':
##
## transpose
##
## The following objects are masked from 'package:lubridate':
##
## hour, isoweek, mday, minute, month, quarter, second, wday, week,
## yday, year
library(ggplot2)
We do not include code for saving figures in this tutorial. If you wish to save figures locally, you can use functions like ggsave() (for ggplot objects) or base R functions such as png(), pdf(), and others.Here are some examples:
## ggsave
# ggsave("my_figure.png", plot = p, width = 8, height = 6, dpi = 300)
## base R function
# png("my_plot.png", width = 800, height = 600, res = 150)
# plot(x, y)
# dev.off()
amerifluxr: https://github.com/chuhousen/amerifluxr# check out data coverage for different sites
full.list = amf_data_coverage(data_product = "BASE-BADM", data_policy = "CCBY4.0")
head(full.list)
## SITE_ID URL
## 1 AR-Bal https://ameriflux.lbl.gov/sites/siteinfo/AR-Bal
## 2 AR-CCa https://ameriflux.lbl.gov/sites/siteinfo/AR-CCa
## 3 AR-CCg https://ameriflux.lbl.gov/sites/siteinfo/AR-CCg
## 4 AR-Cel https://ameriflux.lbl.gov/sites/siteinfo/AR-Cel
## 5 AR-Rom https://ameriflux.lbl.gov/sites/siteinfo/AR-Rom
## 6 AR-TF1 https://ameriflux.lbl.gov/sites/siteinfo/AR-TF1
## publish_years
## 1 2012, 2013
## 2 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020
## 3 2018, 2019, 2020, 2021, 2022, 2023, 2024
## 4
## 5
## 6 2016, 2017, 2018
print(paste("number of sites included in BASE publish is", nrow(full.list)))
## [1] "number of sites included in BASE publish is 761"
# download single site flux/met data
site_name = "US-Syv"
# please keep the codes below commented if you do not have an account
# floc<- amf_download_base(
# user_id = "my_user", # change to your user name
# user_email = "my_email@gmail.com", # change to your email
# site_id = site_name,
# data_product = "BASE-BADM",
# data_policy = "CCBY4.0",
# agree_policy = TRUE,
# intended_use = "other",
# intended_use_text = "Flux-LSM_workshop",
# verbose = TRUE,
# out_dir = tmp_dir
# )
# If you do not have an account, you can directly start from here
tmp_dir <- tempdir(); setwd("~/")
setwd("C:/Users/yl763/Documents/GitHub/FCC_workshop_flux_test")
floc = "AMF_US-Syv_BASE-BADM_29-5.zip"
base <- amf_read_base(file = floc,
unzip = TRUE,
parse_timestamp = TRUE)
# check data coverage for meteorology
# visualizes the BASE data availability for selected AmeriFlux sites, variables, and years.
##PPFD_IN Photosynthetic photon flux density, incoming µmolPhoton m-2 s-1
## TA Air temperature deg C
## RH Relative humidity, range 0-100 %
## VPD Vapor Pressure Deficit hPa
# more variable names: https://ameriflux.lbl.gov/data/aboutdata/data-variables/#base
amf_plot_datayear(site_set = site_name,
var_set = c("PPFD_IN", "TA", "RH", "VPD"),
# year_set = c(2004:2023),
nonfilled_only = TRUE)
# check data coverage for flux data
amf_plot_datayear(site_set = site_name,
var_set = c("FC", "H", "LE"),
# year_set = c(2004:2023),
nonfilled_only = TRUE)
# You can uncomment the line below to check data availability by year for all variables
# amf_plot_datayear(site_set = site_name, nonfilled_only = FALSE)
# subset data for two years
base = base[base$YEAR %in% c(2020,2021), ]
# Recreate other time variables using TIMESTAMP (This is just because I sometimes get error messages in the next step "Intialize EProc" when using the default time variables.)
recreate_time_vars <- function(df) {
df %>%
mutate(
TIMESTAMP = substr(TIMESTAMP_START, 1, 12),
year = substr(TIMESTAMP_START, 1, 4),
month = substr(TIMESTAMP_START, 5, 6),
day = substr(TIMESTAMP_START, 7, 8),
hour = as.numeric(substr(TIMESTAMP_START, 9, 10)) + c(0.5, 1),
date = as.Date(paste(year, month, day, sep = "-")),
doy = yday(date)
)
}
base <- recreate_time_vars(base)
# format the data to be used as input for REddyProc
base_df = data.frame(
'Year' = base$YEAR,
'Hour' = base$hour,
'Date' = base$date,
'Month' = base$month,
'DoY' = base$DOY,
'NEE' = base$FC,
'LE' = base$LE,
'H' = base$LE,
'Rg' = ifelse(base$PPFD_IN_PI_F_1_1_1 < 0, 0, PPFD.to.Rg(base$PPFD_IN_PI_F_1_1_1)),
'Tair' = base$TA_1_1_1,
'Tsoil' = base$TS_1_1_1,
'RH' = ifelse(base$RH_1_1_1 > 100, 100, base$RH_1_1_1),
'VPD' = base$VPD_PI_F_1_1_1,
'Ustar' = base$USTAR_1_1_1,
'TIMESTAMP_START' = as.character(base$TIMESTAMP_START),
'TIMESTAMP_END' = as.character(base$TIMESTAMP_END),
'PPFD' = base$PPFD_IN_PI_F_1_1_1
)
head(base_df)
## Year Hour Date Month DoY NEE LE H Rg Tair
## 1 2020 0.5 2020-01-01 01 1 0.6667798 -4.066935 -4.066935 0 -12.74600
## 2 2020 1.0 2020-01-01 01 1 NA NA NA 0 -13.31267
## 3 2020 1.5 2020-01-01 01 1 NA -2.449585 -2.449585 0 -13.09900
## 4 2020 2.0 2020-01-01 01 1 NA NA NA 0 -12.97467
## 5 2020 2.5 2020-01-01 01 1 0.5290407 -2.824548 -2.824548 0 -13.11433
## 6 2020 3.0 2020-01-01 01 1 NA NA NA 0 -13.14767
## Tsoil RH VPD Ustar TIMESTAMP_START TIMESTAMP_END PPFD
## 1 0.6363515 87.22335 0.2930273 0.2653910 202001010000 202001010030 0
## 2 0.6407161 88.82667 0.2447449 0.2953020 202001010030 202001010100 0
## 3 0.6424356 89.99333 0.2230272 0.3732584 202001010100 202001010130 0
## 4 0.6455239 90.04335 0.2241620 0.6959635 202001010130 202001010200 0
## 5 0.6453573 89.80999 0.2268308 0.8161998 202001010200 202001010230 0
## 6 0.6477322 89.56665 0.2316199 NA 202001010230 202001010300 0
?filterLongRuns
## starting httpd help server ... done
#filterLongRuns : Longer runs, i.e. sequences of numerically identical values, in a series of measurements hint to problems during a noisy measurement, e.g. by sensor malfunction due to freezing. This function, replaces such values in such runs with NA to indicate missing values.
EddyData <- filterLongRuns(base_df, "NEE")
EddyData$Year <- as.numeric(EddyData$Year)
EddyData$Hour <- as.numeric(EddyData$Hour)
EddyData$DoY <- as.numeric(EddyData$DoY)
EddyDataWithPosix <- fConvertTimeToPosix(EddyData, 'YDH', Year = 'Year', Day = 'DoY', Hour = 'Hour')
## Converted time format 'YDH' to POSIX with column name 'DateTime'.
# EProc an R object used in REddyProc with the attributes defined
EProc <- sEddyProc$new(site_name, EddyDataWithPosix, c('NEE','Rg','Tair','VPD', 'Ustar', "H", "LE"))
## New sEddyProc class for site 'US-Syv'
class(EProc)
## [1] "sEddyProc"
## attr(,"package")
## [1] "REddyProc"
# Gapfill meteorological variables using MDS
EProc$sMDSGapFill('Tair', FillAll = FALSE, minNWarnRunLength = NA)
## Initialized variable 'Tair' with 2810 real gaps for gap filling.
## Limited MDS algorithm for gap filling of 'Tair' with LUT(Rg only) and MDC.
## Look up table with window size of 7 days with Rg
## ............................1568
## Mean diurnal course with window size of 0 days: .
## ............0
## Mean diurnal course with window size of 1 days: .
## ............0
## Mean diurnal course with window size of 2 days: .
## ............0
## Look up table with window size of 14 days with Rg
## ............667
## Look up table with window size of 21 days with Rg
## .....575
## Finished gap filling of 'Tair' in 1 seconds. Artificial gaps filled: 35088, real gaps filled: 2810, unfilled (long) gaps: 0.
EProc$sMDSGapFill('Rg', FillAll = FALSE, minNWarnRunLength = NA)
## Initialized variable 'Rg' with 0 real gaps for gap filling.
## Restriced MDS algorithm for gap filling of 'Rg' with no meteo conditions and hence only MDC.
## Finished gap filling of 'Rg' in 0 seconds. Artificial gaps filled: 35088, real gaps filled: 0, unfilled (long) gaps: 0.
EProc$sMDSGapFill('VPD', FillAll = FALSE, minNWarnRunLength = NA)
## Initialized variable 'VPD' with 4270 real gaps for gap filling.
## Limited MDS algorithm for gap filling of 'VPD' with LUT(Rg only) and MDC.
## Look up table with window size of 7 days with Rg
## ..........................................2422
## Mean diurnal course with window size of 0 days: .
## ..................0
## Mean diurnal course with window size of 1 days: .
## ..................0
## Mean diurnal course with window size of 2 days: .
## ..................0
## Look up table with window size of 14 days with Rg
## ..................701
## Look up table with window size of 21 days with Rg
## ...........673
## Look up table with window size of 28 days with Rg
## ....474
## Finished gap filling of 'VPD' in 2 seconds. Artificial gaps filled: 35088, real gaps filled: 4270, unfilled (long) gaps: 0.
# Use MDS to get the "best guess" of true flux
EProc$sMDSGapFill('NEE')
## Initialized variable 'NEE' with 14679 real gaps for gap filling of all 35088 values (to estimate uncertainties).
## Full MDS algorithm for gap filling of 'NEE' with LUT(Rg, VPD, Tair) and MDC.
## Look up table with window size of 7 days with Rg VPD Tair
## ............................................................
## .............................................................
## .............................................................
## .............................................................
## .............................................................
## ...................................................29511
## Look up table with window size of 14 days with Rg VPD Tair
## .......................................................582
## Look up table with window size of 7 days with Rg
## .................................................2982
## Mean diurnal course with window size of 0 days: .
## ....................5
## Mean diurnal course with window size of 1 days: .
## ....................29
## Mean diurnal course with window size of 2 days: .
## ...................10
## Look up table with window size of 21 days with Rg VPD Tair
## ...................1
## Look up table with window size of 28 days with Rg VPD Tair
## ...................0
## Look up table with window size of 35 days with Rg VPD Tair
## ...................0
## Look up table with window size of 42 days with Rg VPD Tair
## ...................0
## Look up table with window size of 49 days with Rg VPD Tair
## ...................0
## Look up table with window size of 56 days with Rg VPD Tair
## ...................0
## Look up table with window size of 63 days with Rg VPD Tair
## ...................0
## Look up table with window size of 70 days with Rg VPD Tair
## ...................0
## Look up table with window size of 14 days with Rg
## ...................770
## Look up table with window size of 21 days with Rg
## ...........685
## Look up table with window size of 28 days with Rg
## .....513
## Finished gap filling of 'NEE' in 29 seconds. Artificial gaps filled: 35088, real gaps filled: 14679, unfilled (long) gaps: 0.
# Calculate residuals and identify outliers
residual <- EProc$sTEMP$NEE_orig - EProc$sTEMP$NEE_fall
IQR <- IQR(residual, na.rm = TRUE)
outlier <- ifelse(abs(residual) > (IQR * 6), 1, 0)
EddieC <- data.frame(
sDateTime = EProc$sTEMP$sDateTime,
NEE_orig = EProc$sTEMP$NEE_orig,
Ustar = EProc$sDATA$Ustar,
NEE_fall = EProc$sTEMP$NEE_fall,
residual = residual,
outlier = outlier
)
# Rename columns
colnames(EddieC) <- c('sDateTime', 'NEE_orig', 'Ustar', 'NEE_fall', 'residual', 'outlier')
# Filter out outliers
EddieC$NEE_filt <- dplyr::if_else(EddieC$outlier > 0, NA_real_, EddieC$NEE_orig)
EddieC$year <- substr(EddieC$sDateTime, 1, 4)
EddieC$doy <- strftime(EddieC$sDateTime, format = "%j")
# Plot the initial outlier detection
EddieC %>%
arrange(as.factor(outlier)) %>%
ggplot(aes(y = NEE_orig, x = as.numeric(doy), color = as.factor(outlier))) +
geom_point(shape = 20, alpha = 0.4) +
theme_minimal() +
labs(x = 'Day of year', y = 'NEE') +
scale_color_manual(values = c('skyblue', 'red')) +
facet_wrap(~year) +
ylim(c(-50, 50)) +
ggtitle("intial outlier detection")
## Warning: Removed 14679 rows containing missing values or values outside the scale range
## (`geom_point()`).
# Re-run the outlier test after initial filtering
EddieC$residual2 <- EddieC$NEE_filt - EddieC$NEE_fall
EddieC$IQR2 <- IQR(EddieC$residual2, na.rm = TRUE)
EddieC$outlier2 <- ifelse(abs(EddieC$residual2) > EddieC$IQR2 * 6, 1, 0)
EddieC$NEE_filt2 <- ifelse(EddieC$outlier2 == 0, EddieC$NEE_filt, NA)
# Plot the re-run outlier detection
EddieC %>%
arrange(as.factor(outlier2)) %>%
ggplot(aes(y = NEE_orig, x = as.numeric(doy), color = as.factor(outlier2))) +
geom_point(shape = 20) +
theme_minimal() +
labs(x = 'Day of year', y = 'NEE') +
scale_color_manual(values = c('skyblue', 'red')) +
facet_wrap(~year) +
ylim(c(-50, 50)) +
ggtitle("re-run outlier detection")
## Warning: Removed 14679 rows containing missing values or values outside the scale range
## (`geom_point()`).
# Remove outliers from the main data
EProc$sDATA$NEE <- EddieC$NEE_filt2
The u* threshold is the minimum u* above which respiration reaches a
plateau (see figure under /plots).
set.seed(2000)
# Here, we are using nSample = 10L for demonstration, please use nSample = 1000L for real research
uStarTh <- EProc$sEstUstarThresholdDistribution(nSample = 10L, probs = c(0.05, 0.5, 0.95))
## Warning in .estimateUStarSeason(...): sEstUstarThreshold: too few finite
## records within season (n = 622) for 7 temperature classes. Need at least n =
## 700. Returning NA for this Season.
## Warning in .estimateUStarSeason(...): sEstUstarThreshold: too few finite
## records within season (n = 633) for 7 temperature classes. Need at least n =
## 700. Returning NA for this Season.
##
## Estimated UStar distribution of:
## uStar 5% 50% 95%
## 1 0.4541172 0.4188225 0.4564736 0.5286447
## by using 10 bootstrap samples and controls:
## taClasses UstarClasses
## 7 20
## swThr minRecordsWithinTemp
## 10 100
## minRecordsWithinSeason minRecordsWithinYear
## 160 3000
## isUsingOneBigSeasonOnFewRecords
## 1
print(uStarTh)
## aggregationMode seasonYear season uStar 5% 50% 95%
## 1 single NA <NA> 0.4541172 0.4188225 0.4564736 0.5286447
## 2 year 2020 <NA> 0.4801723 0.4488349 0.4876796 0.6412940
## 3 year 2021 <NA> 0.4280621 0.3888102 0.4258271 0.4512722
## 4 season 2020 2020001 0.4801723 0.3901247 0.4729320 0.6163569
## 5 season 2020 2020003 0.3993504 0.3877327 0.4484138 0.4839605
## 6 season 2020 2020006 0.2939407 0.2417530 0.3016328 0.3727611
## 7 season 2020 2020009 0.3376365 0.3284783 0.4120114 0.5594615
## 8 season 2021 2020012 0.4280621 0.3950345 0.4302720 0.4641630
## 9 season 2021 2021003 0.4280621 0.3742071 0.4258271 0.4512722
## 10 season 2021 2021006 0.3316664 0.2639192 0.3323809 0.3794716
## 11 season 2021 2021009 0.3632184 0.3288625 0.3656392 0.4183202
## 12 season 2021 2021012 0.4280621 0.3888102 0.4258271 0.4512722
# Define aggregation mode
EProc$sGetUstarScenarios() # by default, annual varying u* is used
## season uStar U05 U50 U95
## 1 2020001 0.4801723 0.4488349 0.4876796 0.6412940
## 2 2020003 0.4801723 0.4488349 0.4876796 0.6412940
## 3 2020006 0.4801723 0.4488349 0.4876796 0.6412940
## 4 2020009 0.4801723 0.4488349 0.4876796 0.6412940
## 5 2020012 0.4280621 0.3888102 0.4258271 0.4512722
## 6 2021003 0.4280621 0.3888102 0.4258271 0.4512722
## 7 2021006 0.4280621 0.3888102 0.4258271 0.4512722
## 8 2021009 0.4280621 0.3888102 0.4258271 0.4512722
## 9 2021012 0.4280621 0.3888102 0.4258271 0.4512722
# Here, we use the single u* threshold
uStar <- uStarTh %>%
dplyr::filter( aggregationMode == "single") %>%
dplyr::select( "uStar", "5%", "50%", "95%")
uStarDf <- cbind(season=na.omit(unique(uStarTh$season)), uStar)
EProc$sSetUstarScenarios(uStarDf)
EProc$sGetUstarScenarios()
## season uStar 5% 50% 95%
## 1 2020001 0.4541172 0.4188225 0.4564736 0.5286447
## 2 2020003 0.4541172 0.4188225 0.4564736 0.5286447
## 3 2020006 0.4541172 0.4188225 0.4564736 0.5286447
## 4 2020009 0.4541172 0.4188225 0.4564736 0.5286447
## 5 2020012 0.4541172 0.4188225 0.4564736 0.5286447
## 6 2021003 0.4541172 0.4188225 0.4564736 0.5286447
## 7 2021006 0.4541172 0.4188225 0.4564736 0.5286447
## 8 2021009 0.4541172 0.4188225 0.4564736 0.5286447
## 9 2021012 0.4541172 0.4188225 0.4564736 0.5286447
# if you want to use seasonal or annual varying u*:
# EProc$useSeaonsalUStarThresholds()# seasonal varying u*
# EProc$useAnnualUStarThresholds() # annual varying u*
# EProc$sPlotNEEVersusUStarForSeason(format = "pdf") # keep this line commented when use binder, then check the figure under /plots on github repo
Bonus training: create a figure comparing u* threshold determined using different aggregationMode.
Fluxes are expected similar if they are:
at similar environmental conditions
and close in time
Fluxes are expected similar if they are at the same time of the day
and not too many days away;
The methods are effective for short gaps as the missing values are
replaced by the average of response variables under similar weather
conditions in a small-time window;
Quality flag increases with fewer variables and larger time windows:
# fingerplot: inspect gaps needed to be filled
EProc$sPlotFingerprintY("NEE", Year = 2020); EProc$sPlotFingerprintY("NEE", Year = 2021)
# use MDS to gapfill flux data
EProc$sMDSGapFillUStarScens('NEE')
## Ustar filtering (u * Th_1 = 0.45411720625), marked 33% of the data as gap
## Initialized variable 'NEE' with 20978 real gaps for gap filling of all 35088 values (to estimate uncertainties).
## Full MDS algorithm for gap filling of 'NEE.Ustar_uStar_fqc_0' with LUT(Rg, VPD, Tair) and MDC.
## Look up table with window size of 7 days with Rg VPD Tair
## ............................................................
## .............................................................
## .............................................................
## .............................................................
## .............................................................
## ...................................................29068
## Look up table with window size of 14 days with Rg VPD Tair
## ............................................................
## .905
## Look up table with window size of 7 days with Rg
## ...................................................3092
## Mean diurnal course with window size of 0 days: .
## ....................5
## Mean diurnal course with window size of 1 days: .
## ....................29
## Mean diurnal course with window size of 2 days: .
## ...................10
## Look up table with window size of 21 days with Rg VPD Tair
## ...................1
## Look up table with window size of 28 days with Rg VPD Tair
## ...................0
## Look up table with window size of 35 days with Rg VPD Tair
## ...................0
## Look up table with window size of 42 days with Rg VPD Tair
## ...................0
## Look up table with window size of 49 days with Rg VPD Tair
## ...................0
## Look up table with window size of 56 days with Rg VPD Tair
## ...................0
## Look up table with window size of 63 days with Rg VPD Tair
## ...................0
## Look up table with window size of 70 days with Rg VPD Tair
## ...................0
## Look up table with window size of 14 days with Rg
## ...................768
## Look up table with window size of 21 days with Rg
## ............681
## Look up table with window size of 28 days with Rg
## .....529
## Finished gap filling of 'NEE' in 28 seconds. Artificial gaps filled: 35088, real gaps filled: 20978, unfilled (long) gaps: 0.
## Ustar filtering (u * Th_1 = 0.418822541913729), marked 31% of the data as gap
## Initialized variable 'NEE' with 20411 real gaps for gap filling of all 35088 values (to estimate uncertainties).
## Full MDS algorithm for gap filling of 'NEE.Ustar_5%_fqc_0' with LUT(Rg, VPD, Tair) and MDC.
## Look up table with window size of 7 days with Rg VPD Tair
## ............................................................
## .............................................................
## .............................................................
## .............................................................
## .............................................................
## ...................................................29158
## Look up table with window size of 14 days with Rg VPD Tair
## ...........................................................843
## Look up table with window size of 7 days with Rg
## ..................................................3064
## Mean diurnal course with window size of 0 days: .
## ....................5
## Mean diurnal course with window size of 1 days: .
## ....................29
## Mean diurnal course with window size of 2 days: .
## ...................10
## Look up table with window size of 21 days with Rg VPD Tair
## ...................1
## Look up table with window size of 28 days with Rg VPD Tair
## ...................0
## Look up table with window size of 35 days with Rg VPD Tair
## ...................0
## Look up table with window size of 42 days with Rg VPD Tair
## ...................0
## Look up table with window size of 49 days with Rg VPD Tair
## ...................0
## Look up table with window size of 56 days with Rg VPD Tair
## ...................0
## Look up table with window size of 63 days with Rg VPD Tair
## ...................0
## Look up table with window size of 70 days with Rg VPD Tair
## ...................0
## Look up table with window size of 14 days with Rg
## ...................768
## Look up table with window size of 21 days with Rg
## ............681
## Look up table with window size of 28 days with Rg
## .....529
## Finished gap filling of 'NEE' in 29 seconds. Artificial gaps filled: 35088, real gaps filled: 20411, unfilled (long) gaps: 0.
## Ustar filtering (u * Th_1 = 0.456473581875), marked 33% of the data as gap
## Initialized variable 'NEE' with 21021 real gaps for gap filling of all 35088 values (to estimate uncertainties).
## Full MDS algorithm for gap filling of 'NEE.Ustar_50%_fqc_0' with LUT(Rg, VPD, Tair) and MDC.
## Look up table with window size of 7 days with Rg VPD Tair
## ............................................................
## .............................................................
## .............................................................
## .............................................................
## .............................................................
## ...................................................29062
## Look up table with window size of 14 days with Rg VPD Tair
## ............................................................
## .909
## Look up table with window size of 7 days with Rg
## ...................................................3094
## Mean diurnal course with window size of 0 days: .
## ....................5
## Mean diurnal course with window size of 1 days: .
## ....................29
## Mean diurnal course with window size of 2 days: .
## ...................10
## Look up table with window size of 21 days with Rg VPD Tair
## ...................1
## Look up table with window size of 28 days with Rg VPD Tair
## ...................0
## Look up table with window size of 35 days with Rg VPD Tair
## ...................0
## Look up table with window size of 42 days with Rg VPD Tair
## ...................0
## Look up table with window size of 49 days with Rg VPD Tair
## ...................0
## Look up table with window size of 56 days with Rg VPD Tair
## ...................0
## Look up table with window size of 63 days with Rg VPD Tair
## ...................0
## Look up table with window size of 70 days with Rg VPD Tair
## ...................0
## Look up table with window size of 14 days with Rg
## ...................768
## Look up table with window size of 21 days with Rg
## ............681
## Look up table with window size of 28 days with Rg
## .....529
## Finished gap filling of 'NEE' in 28 seconds. Artificial gaps filled: 35088, real gaps filled: 21021, unfilled (long) gaps: 0.
## Ustar filtering (u * Th_1 = 0.528644729070933), marked 37% of the data as gap
## Initialized variable 'NEE' with 22024 real gaps for gap filling of all 35088 values (to estimate uncertainties).
## Full MDS algorithm for gap filling of 'NEE.Ustar_95%_fqc_0' with LUT(Rg, VPD, Tair) and MDC.
## Look up table with window size of 7 days with Rg VPD Tair
## ............................................................
## .............................................................
## .............................................................
## .............................................................
## .............................................................
## ...................................................28810
## Look up table with window size of 14 days with Rg VPD Tair
## ............................................................
## ...1076
## Look up table with window size of 7 days with Rg
## ....................................................3178
## Mean diurnal course with window size of 0 days: .
## ....................5
## Mean diurnal course with window size of 1 days: .
## ....................29
## Mean diurnal course with window size of 2 days: .
## ...................10
## Look up table with window size of 21 days with Rg VPD Tair
## ...................1
## Look up table with window size of 28 days with Rg VPD Tair
## ...................0
## Look up table with window size of 35 days with Rg VPD Tair
## ...................0
## Look up table with window size of 42 days with Rg VPD Tair
## ...................0
## Look up table with window size of 49 days with Rg VPD Tair
## ...................0
## Look up table with window size of 56 days with Rg VPD Tair
## ...................0
## Look up table with window size of 63 days with Rg VPD Tair
## ...................0
## Look up table with window size of 70 days with Rg VPD Tair
## ...................0
## Look up table with window size of 14 days with Rg
## ...................768
## Look up table with window size of 21 days with Rg
## ............681
## Look up table with window size of 28 days with Rg
## .....530
## Finished gap filling of 'NEE' in 28 seconds. Artificial gaps filled: 35088, real gaps filled: 22024, unfilled (long) gaps: 0.
# gaps in flux data are filled
EProc$sPlotFingerprintY("NEE_50._f", Year = 2020); EProc$sPlotFingerprintY("NEE_50._f", Year = 2021)
Reco is modeled and GPP is computed as Reco - NEE;
There are three approaches to partition FCO₂ available in REddyProc:
Night-time partitioning (NT): This approach relies on the assumption that fluxes measured at night consist purely of respiration. It uses nighttime fluxes to estimate a seasonally varying reference respiration rate and its sensitivity to temperature (Reichstein 2005).
Day-time partitioning (DT): This method uses a light response curve based on daytime data to estimate the reference respiration rate, and then applies a temperature response function to estimate Reco (Lasslop 2010).
Modified day-time partitioning: Both NT and DT assume that the same Rref is applicable during daytime as at night—a questionable assumption due to the potential for inhibition of leaf respiration in the light (Keenan 2019).
# TimeZoneHour: time zone offset from UTC (Coordinated Universal Time), without daytime saving. This is a site-specific input.
EProc$sSetLocationInfo(LatDeg = 46.2420, LongDeg =-89.3477, TimeZoneHour = -6)
EProc$sFillVPDFromDew() # fill longer gaps still present in VPD_f
EProc$sMRFluxPartitionUStarScens("NEE") # night-time
## Start flux partitioning for variable NEE with temperature Tair_f.
## Estimate of the temperature sensitivity E_0 from short term data: 84.51.
## Regression of reference temperature R_ref for 169 periods.
## Start flux partitioning for variable NEE with temperature Tair_f.
## Estimate of the temperature sensitivity E_0 from short term data: 84.51.
## Regression of reference temperature R_ref for 169 periods.
## Start flux partitioning for variable NEE with temperature Tair_f.
## Estimate of the temperature sensitivity E_0 from short term data: 84.51.
## Regression of reference temperature R_ref for 169 periods.
## Start flux partitioning for variable NEE with temperature Tair_f.
## Estimate of the temperature sensitivity E_0 from short term data: 84.51.
## Regression of reference temperature R_ref for 169 periods.
# If you want to expolore day-time or night-time partitioning
# EProc$sGLFluxPartitionUStarScens() # day-time
# EProc$sTKFluxPartitionUStarScens() # modified day-time
names(EProc$sTEMP)
## [1] "sDateTime" "Tair_orig" "Tair_f"
## [4] "Tair_fqc" "Tair_fall" "Tair_fall_qc"
## [7] "Tair_fnum" "Tair_fsd" "Tair_fmeth"
## [10] "Tair_fwin" "Rg_orig" "Rg_f"
## [13] "Rg_fqc" "Rg_fall" "Rg_fall_qc"
## [16] "Rg_fnum" "Rg_fsd" "Rg_fmeth"
## [19] "Rg_fwin" "VPD_orig" "VPD_f"
## [22] "VPD_fqc" "VPD_fall" "VPD_fall_qc"
## [25] "VPD_fnum" "VPD_fsd" "VPD_fmeth"
## [28] "VPD_fwin" "NEE_orig" "NEE_f"
## [31] "NEE_fqc" "NEE_fall" "NEE_fall_qc"
## [34] "NEE_fnum" "NEE_fsd" "NEE_fmeth"
## [37] "NEE_fwin" "season" "Ustar_uStar_Thres"
## [40] "Ustar_uStar_fqc" "NEE_uStar_orig" "NEE_uStar_f"
## [43] "NEE_uStar_fqc" "NEE_uStar_fall" "NEE_uStar_fall_qc"
## [46] "NEE_uStar_fnum" "NEE_uStar_fsd" "NEE_uStar_fmeth"
## [49] "NEE_uStar_fwin" "Ustar_5._Thres" "Ustar_5._fqc"
## [52] "NEE_5._orig" "NEE_5._f" "NEE_5._fqc"
## [55] "NEE_5._fall" "NEE_5._fall_qc" "NEE_5._fnum"
## [58] "NEE_5._fsd" "NEE_5._fmeth" "NEE_5._fwin"
## [61] "Ustar_50._Thres" "Ustar_50._fqc" "NEE_50._orig"
## [64] "NEE_50._f" "NEE_50._fqc" "NEE_50._fall"
## [67] "NEE_50._fall_qc" "NEE_50._fnum" "NEE_50._fsd"
## [70] "NEE_50._fmeth" "NEE_50._fwin" "Ustar_95._Thres"
## [73] "Ustar_95._fqc" "NEE_95%_orig" "NEE_95%_f"
## [76] "NEE_95%_fqc" "NEE_95%_fall" "NEE_95%_fall_qc"
## [79] "NEE_95%_fnum" "NEE_95%_fsd" "NEE_95%_fmeth"
## [82] "NEE_95%_fwin" "VPDfromDew" "PotRad_5%"
## [85] "FP_NEEnight_5%" "FP_Temp_5%" "E_0_5%"
## [88] "R_ref_5%" "Reco_5%" "GPP_5%_f"
## [91] "GPP_5%_fqc" "PotRad_50%" "FP_NEEnight_50%"
## [94] "FP_Temp_50%" "E_0_50%" "R_ref_50%"
## [97] "Reco_50%" "GPP_50%_f" "GPP_50%_fqc"
## [100] "PotRad_95%" "FP_NEEnight_95%" "FP_Temp_95%"
## [103] "E_0_95%" "R_ref_95%" "Reco_95%"
## [106] "GPP_95%_f" "GPP_95%_fqc" "PotRad_uStar"
## [109] "FP_NEEnight_uStar" "FP_Temp_uStar" "E_0_uStar"
## [112] "R_ref_uStar" "Reco_uStar" "GPP_uStar_f"
## [115] "GPP_uStar_fqc"
Bonus training: comparing different partitioning method;
FilledEddyData <- EProc$sExportResults()
combined.df <- cbind(EddyData, FilledEddyData)
names(combined.df)
## [1] "Year" "Hour" "Date"
## [4] "Month" "DoY" "NEE"
## [7] "LE" "H" "Rg"
## [10] "Tair" "Tsoil" "RH"
## [13] "VPD" "Ustar" "TIMESTAMP_START"
## [16] "TIMESTAMP_END" "PPFD" "Tair_orig"
## [19] "Tair_f" "Tair_fqc" "Tair_fall"
## [22] "Tair_fall_qc" "Tair_fnum" "Tair_fsd"
## [25] "Tair_fmeth" "Tair_fwin" "Rg_orig"
## [28] "Rg_f" "Rg_fqc" "Rg_fall"
## [31] "Rg_fall_qc" "Rg_fnum" "Rg_fsd"
## [34] "Rg_fmeth" "Rg_fwin" "VPD_orig"
## [37] "VPD_f" "VPD_fqc" "VPD_fall"
## [40] "VPD_fall_qc" "VPD_fnum" "VPD_fsd"
## [43] "VPD_fmeth" "VPD_fwin" "NEE_orig"
## [46] "NEE_f" "NEE_fqc" "NEE_fall"
## [49] "NEE_fall_qc" "NEE_fnum" "NEE_fsd"
## [52] "NEE_fmeth" "NEE_fwin" "season"
## [55] "Ustar_uStar_Thres" "Ustar_uStar_fqc" "NEE_uStar_orig"
## [58] "NEE_uStar_f" "NEE_uStar_fqc" "NEE_uStar_fall"
## [61] "NEE_uStar_fall_qc" "NEE_uStar_fnum" "NEE_uStar_fsd"
## [64] "NEE_uStar_fmeth" "NEE_uStar_fwin" "Ustar_5._Thres"
## [67] "Ustar_5._fqc" "NEE_5._orig" "NEE_5._f"
## [70] "NEE_5._fqc" "NEE_5._fall" "NEE_5._fall_qc"
## [73] "NEE_5._fnum" "NEE_5._fsd" "NEE_5._fmeth"
## [76] "NEE_5._fwin" "Ustar_50._Thres" "Ustar_50._fqc"
## [79] "NEE_50._orig" "NEE_50._f" "NEE_50._fqc"
## [82] "NEE_50._fall" "NEE_50._fall_qc" "NEE_50._fnum"
## [85] "NEE_50._fsd" "NEE_50._fmeth" "NEE_50._fwin"
## [88] "Ustar_95._Thres" "Ustar_95._fqc" "NEE_95%_orig"
## [91] "NEE_95%_f" "NEE_95%_fqc" "NEE_95%_fall"
## [94] "NEE_95%_fall_qc" "NEE_95%_fnum" "NEE_95%_fsd"
## [97] "NEE_95%_fmeth" "NEE_95%_fwin" "VPDfromDew"
## [100] "PotRad_5%" "FP_NEEnight_5%" "FP_Temp_5%"
## [103] "E_0_5%" "R_ref_5%" "Reco_5%"
## [106] "GPP_5%_f" "GPP_5%_fqc" "PotRad_50%"
## [109] "FP_NEEnight_50%" "FP_Temp_50%" "E_0_50%"
## [112] "R_ref_50%" "Reco_50%" "GPP_50%_f"
## [115] "GPP_50%_fqc" "PotRad_95%" "FP_NEEnight_95%"
## [118] "FP_Temp_95%" "E_0_95%" "R_ref_95%"
## [121] "Reco_95%" "GPP_95%_f" "GPP_95%_fqc"
## [124] "PotRad_uStar" "FP_NEEnight_uStar" "FP_Temp_uStar"
## [127] "E_0_uStar" "R_ref_uStar" "Reco_uStar"
## [130] "GPP_uStar_f" "GPP_uStar_fqc"
# Have a look at the gapfilled data
# diurnal pattern of FC
ggplot(combined.df, aes(x = Hour, y = NEE_50._f)) +
geom_point(col = "grey") +
geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
stat_summary(fun = mean, geom = "line") +
facet_wrap(~ interaction(Month, Year)) +
theme_minimal() +
ggtitle("FC") + ylab(expression(FCO[2] ~ "(" * mu * "mol" ~ m^{-2} ~ s^{-1} * ")"))
# diurnal pattern of GPP
ggplot(combined.df, aes(x = Hour, y = `GPP_50%_f`)) +
geom_point(col = "grey") +
stat_summary(fun = mean, geom = "line") +
geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
facet_wrap(~ interaction(Month, Year)) +
theme_minimal() +
ggtitle("GPP") + ylab(expression(GPP ~ "(" * mu * "mol" ~ m^{-2} ~ s^{-1} * ")"))
## Warning: Removed 14906 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 14906 rows containing missing values or values outside the scale range
## (`geom_point()`).
# diural pattern of Reco
ggplot(combined.df, aes(x = Hour, y = `Reco_50%`)) +
geom_point(col = "grey") +
stat_summary(fun = mean, geom = "line") +
geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
facet_wrap(~ interaction(Month, Year)) +
theme_minimal() +
ggtitle("Reco") + ylab(expression(Reco ~ "(" * mu * "mol" ~ m^{-2} ~ s^{-1} * ")"))
Bonus training: gap-fill H or LE, and plot the diurnal patterns.